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A cluster algorithm formulated in continuous (imaginary) time is presented for Ising models 
in a transverse field. It works directly with an infinite number of time-slices in the imaginary 
time direction, avoiding the necessity to take this limit explicitly. The algorithm is tested at the 
zero-temperature critical point of the pure two-dimensional (2d) transverse Ising model. Then it is 
applied to the 2d Ising ferromagnet with random bonds and transverse fields, for which the phase 
diagram is determined. Finite size scaling at the quantum critical point as well as the study of the 
quantum Griffiths-McCoy phase indicate that the dynamical critical exponent is infinite as in Id. 



Quantum phase transitions in random transverse Ising 
models at and close to their quantum critical point at 
zero temperature (T = 0) have attracted a lot of in- 
terest recentlya. In particular in one dimension many 
astonishing results have been obtained with powerful an- 
alytical (real space renormalization groupn, field theorytl) 
and numerical (exact diagonalization after mapping on a 
free fermion modelotj) tools. Among the most important 
results is that at the critical point time scales diverge ex- 
ponentially fast, implying that the dynamical exponent is 
z cr j t = oo. Such a behavior is reminiscent of thermally ac- 
tivated dynamics in classical random field systemsB, and 
has also been proposed to be present at the Anderson- 
Mott transition of disordered electronic systems at the 
quantum phase transition!! 

The other important result for random transverse Ising 
models in one dimension is that close to the quantum crit- 
ical point there is a whole region in which various suscep- 
tibilities diverge for T — > and that all thcsfi|-singular- 
ities, also called Griffiths-McCoy singularitiesaiJ, can be 
parameterized by a single dynamical exponent z(S) that 
varies continuously with the difference 5 to the critical 
pont and diverges for S — > 0. It is not clear in how far 
these properties, i.e. z cr ;t = oo and z(S — * 0) = z cr it also 
apply in higher dimensions: In 2 and 3-dimcnsional trans- 
verse Ising spin, glass models a finite value for z cr it has 
been reported£3 and for finite .dimensional bond diluted 
ferromagnets it has been shownliil that z CT n is infinite only 
at the percolation threshold. 

In this letter we therefore consider the two-dimensional 
random ferromagnet (without dilution) in a transverse 
field with the help of a new Monte-Carlo cluster algo- 
rithm that is particularly suited to handle the inher- 
ent difficulties in the study of such a random quantum 
system: The origin of the Griffiths-McCoy singularities 
are strongly coupled clusters (with strong ferromagnetic 
bonds and weak transverse fields) which are extremely 



hard to equilibrate in a conventional Monte-Carlo algo- 
rithm, therefore we had to use a cluster method. More- 
over, the exponent z(5) is a non-universal quantity for 
which reason we really have to perform the so-called 
Trotter-limiL explained below. In a related work, Pich 
and Youngtj have studied essentially the same model 
with a discrete time cluster algorithm, not performing 
the Trotter limit and therefore concentrating on the crit- 
ical behavior. 

A continuous time algorithm that incorporates—thus 
limit right from the beginning (in the spirit of ref.EirEj) 
is the most efficient method we can think of. In the first 
part of this letter we present the method we use and ap- 
ply it to the pure case in two dimensions, in the second 
part we present our results for the random case. Details 
and additional results will be published elsewheret3. 

The system we are interested in is defined by the quan- 
tum mechanical Hamiltonian 



(i) 



in which <Xj are spin-i operators located on any d- 
dimensional lattice (belw we consider aLxL square lat- 
tice with periodic boundary conditions), (ij) indicate all 
nearest neighbor pairs on this particular lattice, Jij > 
ferromagnetic interactions (either uniform Jy = J or 
random) and Tj transverse fields (either uniform I\ = T 
or random - note that the sign of T% can always be gauged 
away by a local spin rotation) . 

To derive our continuous time algorithm we use the Su- 
zuki-Trotter decomposition to represent the free energy 
T of the system ([!]) at inverse temperature (3 — l/T^as 
the limit of a (d+l)-dimensional classical Ising modeO: 
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Here the additional index r = 1, . . . , L T of the now clas- 
sical Ising spin variables Si(r) = ±1 labels the L T d- 
dimensional (imaginary) time slices within which spins 



interact via K, 



At J,-. 



Jij and among which they inter- 
act with strength K\ = — i In tanh Arr^ . The number of 
time slices L T is related to Ar by At = (3/L T , so that 
the limit At — > in (Q) implies L T — > oo. 

Taking the limit At — > consecutive spins with 
the same value along the imaginary time direction, e.g. 
Si(r)=Si(T + 1)=- • -=Si(T + N), form continuous seg- 
ments Si{[r,T + t]} of length t — N ■ At rather than 
individual lattice points. Since we are going to take this 
limit implicitly we will have to consider these imaginary 
time segments as the dynamical objects in a Monte-Carlo 
algorithm, and not the individual spin values at discrete 
imaginary times any more. These segments correspond 
to continuous, uninterrupted pieces of a spin's world line 
during which it has the same value, say +1, and its two 
ends, in the following called cuts, are those times when 
this spin changes to another value, say — 1. 

As the next step we apply the-scheme of the Swendsen- 
Wang cluster update methodEJ within the aforemen- 
tioned implicit continuous time limit. Remember that 
in this method in order to construct the clusters to be 
flipped at random, neighboring spins pointing in the 
same direction are connected with a certain probabil- 
ity p: for neighbors in the space direction, for instance 
Si(r) and Sj(r) with (ij) being nearest neighbors, it is 
Pij = 1 — exp(— IKij) = 2ArJij +C(At 2 ) and for neigh- 
bors in imaginary time, for instance Si(r) and Si(r+1), it 
is p\ = 1 - exp(-2Kl) = 1 - AtT, + C(At 2 ). These spin 
connection probabilities are now translated into proba- 
bilities for creating (cutting) and connecting segments. 

The probability for connecting spins along the imagi- 
nary time direction at a particular site i over a finite time 
interval of length t < (3 is given by the probability to set 
tj At bonds, i.e. in the limit At — > 
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This means for each site one generates new cuts in ad- 
dition to the old ones from the already existing seg- 
ments via a Poisson process with decay time 1 /Ti along 
the imaginary time direction. Next one connects seg- 
ments on neighboring sites i and j that have the same 
state and a nonvanishing time-overlap t, e.g. Si{[ti,t2]} 
and Sj { [£3 , ti\ } with t now the length of the interval 
[£1,^2] H [£3, £4]). The probability for not connecting the 
two segments is given by 



(1-Py) t/Ar = (1-2AtJ«) 



t/AT 



exp(-2Jyi) (4) 



in the limit At — > 0. Finally, one identifies clusters of 
connected segments and assigns each of them (and here- 
with all of the segments belonging to it) one value +1 
or —1 with equal probability. In Fig. [j] the individ- 
ual steps of this cluster update procedure are illustrated. 
The measurement of observables is straightforward: for 



instance the local magnetization at site i for one par- 
ticular configuration of segments is simply given by the 
difference between the total length of all +segments and 
the total length of all —segments, divided by j3. The 
expectation value for the local susceptibility is 
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where (• • -)qm is the quantum mechanical expectation 
value for model ([!]) and (• • -}mc is an average over all 
configurations generated during the Monte-Carlo run. 
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FIG. 1. Sketch of the cluster construction procedure: a) 
Configuration of segments before update — full (broken) lines 
correspond to worldline segments with spin up (down), b) In- 
sertion of new cuts (crosses) in addition to old ones (bars) 
according to a Poisson process along each worldline described 
by ©• c ) Connection of segments with probabilities given by 
(M). d) Configuration of segments after assigning randomly a 
spin value to each cluster of connected segments. Before start- 
ing a new cycle (with A) the redundant cuts within segments 
are removed. 

For an actual implementation of the algorithm one 
has to provide sufficient memory in which the informa- 
tion about the segments Si{[t\, £2]} (site index i, starting 
point t\, end point t 2 and state +1 or —1) is stored in 
linked lists. The number of segments, which of course 
fluctuates, is typically of the order 0{(3T max L d ), whereas 
in discrete time algorithm the number of spin values to 
be stored is (3L d /Ar, which diverges in the limit At — > 0. 

We tested our code first for the one-dimensional pure 
case (d = 1, = 1 and Ti = 7), for which the critical 
properties at T = are identical to those of the classical 
2-dimcnsional Ising model and are well known exactlycSl 

(r c = i,i/ = i,/? = i/8,z = i). 

As a first non-trivial application we considered the pure 
case (again = 1 and = T) in two dimensions, where 
the quantum critical point at T c is in the universality 
class of the classical 3-dimensional Ising model. Close 
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to this quantum critical point quantities are expected to 
obey finite size scaling forms 

(<D)=L x °q(L 1 /»5,L*/f3) (6) 

with 5 = r — T c and xq the finite size scaling expo- 
nent of the quantity O. From the equivalence with the 
3-dimensional classical Ising model one knows the dy- 



namical exponent a priori to be z = 1, thus we can per- 
form conventional one-parameter finite size scaling if one 
chooses the aspect ratio (3/L to be constant. In Tab. I 
we list the results of the finite size scaling analysis of our 
data for system sizes up to L = 32 (and (3 = 8). Our 
estimates for the critical field value T c = 3.044(1) and 
the thermal exponent^ = 0.625(5) agree well with— the 
series expansion resulto and a recent DMRG studyEl. 
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TABLE I. Estimates for the critical field strength and the critical exponents for the quantum phase transition of the pure 
two-dimensional transverse Ising model. The left column indicates the quantity for which the finite size scaling analysis with 
constant aspect ratio /3 = L/4 has been performed (m = magnetization, \ — uniform susceptibility and g = dimensionless ratio 
of moments). The values in brackets are obtained from the scaling relation 7/1/ = d + z — 2(3 jv id — 2 and z = 1 here). 



The real challenge for our algorithm (and our main mo- 
tivation for implementing it) is the random case, which 
we consider now. The ferromagnetic couplings as well as 
the transverse field are now quenched random variables, 
which wc define both as being uniformly distributed, i.e. 



P(Jij) = 
P'Fi) = 



1, for < Jij < 1 
0, otherwise 

I?" 1 , for < Ti < r 
0, otherwise 



(7) 



Observables now have to be averaged over a large num- 
ber (10 3 -10 4 ) of disorder configurations, which will be 
indicated by [. . .] av . 




FIG. 2. The phase diagram of the two-dimensional ran- 
dom transverse Ising model. PM means paramagnetic, FM 
means ferromagnetic, T c = 1.05(1) is the critical tempera- 
ture of the classical random Ising ferromagnet with a uniform 
bond distribution, and P c = 4.2(2) the location of the quan- 
tum critical point we are interested in. 



The main interest in our investigation is the determi- 
nation of the dynamical exponent at and above the quan- 
tum critical point. Since both are unknown, it is hard to 
work with a finite size scaling form like (0) . Therefore we 
have chosen the following procedure: first we determine 
the (T,r) phase diagram of the model (Q) with (Q) by 
calculating the averaged ratio of moments [g] av for var- 
ious system sizes L at fixed temperature T > and T. 
At any finite temperature the system is expected to be 
in the universality class of the classical two-dimensional 
classical random bond Ising ferromagnet. For small tem- 
peratures the crossover region to the quantum universal- 
ity class becomes larger, necessitating large system sizes 
(we went up to L — 32) to get a reliable estimate of 
r c (T). By extrapolating the latter to T = 0, see Fig. ||, 
we obtain for the location of the quantum critical point 
r c = r c (T = 0) = 4.2 ±0.2. 

For r = T c it is 5 = and any observable is ex- 
pected to be a function of the aspect ratio P/L z alone. 
On the other hand, if the two-dimensional case, which 
we consider here, is similar to the one-dimensional c 
where one has unconventional (or activated) scaling! 
one would expect z = 00, which implies that ln/3/L^ 
would be the appropriate scaling variable (with ip a fit- 
parameter, playing the rolo-,af-the barrier exponent in 
classical activated dynamicsQl3'E2l). Our data scale much 
better in according to the latter scenaria with ip close to 
0.5 as in Id (details will be published cordance 
with the results reported by Pich and Youngcj. 

Now we turn our attention to the Griffiths-McCoy re- 
gion in the disordered phase (r > T c ). Due to the pres- 
ence of strongly coupled regions in the system the prob- 
ability distribution of excitation energies (essentially in- 
verse tunneling times for these ferromagnetically ordered 
clusters) becomes extremely broad. As a consequence we 
expect the probability distribution of local susceptibili- 
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ties to have an algebraic tail at T = 




^ (In Xiocal ) ~ - -^py In Xiocal (8) 

where f2(lnxi oca i) is the probability for the logarithm 
of the local susceptibility Xi at site i to be larger than 
In xiocal- The dynamical exponent z(T) varies continu- 
ously with the distance from the critical point and pa- 
rameterizes the strengths of the Griffiths-McCoy singu- 
larities also present in other observables. At finite tem- 
peratures f2 is chopped off at j3, and close to the critical 
point one expects finite size corrections as long as L or /3 
are smaller than the spatial correlation length or imag- 
inary correlation time, respectively. We used (3 < 1000 
and averaged over at least 512 samples. 



Here we found indications for an exponential divergence 
of time scales (z = oo) and also the dynamical exponent 
parameterizing the strength of the quantum Griffiths- 
McCoy singularities extrapolates to z(T — > r c ) = oo. In 
this respect the phenomenology of the one-dimensional 
model seems to extend to higher dimension. Another 
aspect is the different scaling behavior of average and 
typical correlations in tfijC-|One-dimensional case, which 
seems to be present in 2dli3, too. 
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Forschungsgemeinschaft (DFG) and is grateful to the 
Toho University Department of Physics for kind hospi- 
tality. N.K.'s work was supported by the grant-in-aid 
(No. 09740320) from the ministry of education, science 
and culture. 
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FIG. 3. The value of d/z(T) obtained from analyzing the 
integrated probability distribution of In Xiocal according to (^) 
in the Griffiths-McCoy region. The vertical line indicates the 
(approximate) region of the critical point at r c ~ 4.2, the 
open circle corresponds to z(T c ) = oo and the horizontal line 
at d/z — 1 indicates the expected limit limr^oo z(T) = d. 
The broken line is just a guide to the eye. 

In Fig. H we show our result for d/z(T) in the Griffiths- 
McCoy region. For T — > oo we expect d/z(T) — » 1, 
since this is the results for isolated spins in random 
fields with non- vanishing probability weight at Tj = 0. 
The more interesting limit is T — * T c . Due to the 
aforementioned finite size effects an accurate estimate of 
d/z becomes increasingly hard approaching the critical 
point. Nevertheless, the data are well compatible with 
limr^r c Z (T) = oo, implying that here, in analogy to the 
one-dimensional caseaij this limit and the critical dynam- 
ical exponent z cr ;t agree. 

To summarize we have presented a new Monte-Carlo 
cluster algorithm in continuous imaginary time with 
which we studied the random transverse Ising model 
in two dimensions. We determined the temperature- 
transverse field phase diagram, estimated the location of 
the quantum critical point at zero temperature and per- 
formed a finite size scaling analysis at the critical point. 
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